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Abstract. We introduce a family of fourth order two-step methods that preserve the energy 
function of canonical polynomial Hamiltonian systems. Each method in the family may be viewed 
as a correction of a linear two-step method, where the correction term is 0(/i 5 ) (h is the stepsize of 
integration). The key tools the new methods are based upon are the line integral associated with 
a conservative vector field (such as the one denned by a Hamiltonian dynamical system) and its 
discretization obtained by the aid of a quadrature formula. Energy conservation is equivalent to the 
requirement that the quadrature is exact, which turns out to be always the case in the event that 
the Hamiltonian function is a polynomial and the degree of precision of the quadrature formula is 
high enough. The non-polynomial case is also discussed and a number of test problems are finally 
presented in order to compare the behavior of the new methods to the theoretical results. 
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1. Introduction and Background. We consider canonical Hamiltonian sys- 
tems in the form 

^ = JVH(y), J=(_° Im 7 o n )> f(*o) = Va G M 2m , (1-1) 

where H(y) is a smooth real- valued function. Our interest is in researching numerical 
methods that provide approximations y n ~ y(to+nh) to the true solution along which 
the energy is precisely conserved, namely 

H(y n ) = H(yo), for all stepsizes h < ho. (1-2) 

The study of energy-preserving methods form a branch of geometrical numeri- 
cal integration, a research topic whose main aim is preserving qualitative features of 
simulated differential equations. In this context, symplectic methods have had con- 
siderable attention due to their good long-time behavior as compared to standard 
methods for ODEs [25, 15, 21]. A related interesting approach based upon exponen- 
tial/trigonometric fitting may be found in [20, 27, 26]. Unfortunately, symplecticity 
cannot be fully combined with the energy preservation property [16], and this partly 
explains why the latter has been absent from the scene for a long time. 

Among the first examples of energy-preserving methods we mention discrete gra- 
dient schemes [17, 23] which are defined by devising discrete analogs of the gradient 
function. The first formulae in this class had order at most two but recently discrete 
gradient methods of arbitrarily high order have been researched by considering the 
simpler case of systems with one-degree of freedom [12, 13]. 

Here, the key tool we wish to exploit is the well-known line integral associated 
with conservative vector fields, such us the one defined at (1.1), as well as its discrete 
version, the so called discrete line integral. Interestingly, the line integral provides a 
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means to check the energy conservation property, namely 



H{y{tx)) - H{y ) =f VH(y)dy = h f y 1 (t + rh) T X7 H{y{t + Th))Ar 
Jyo->v(ti) JO 



= h f V T H(y(t + Th))J T VH{y{t a + r/i))dr = 0, 
Jo 

with h = t\ — to, that can be easily converted into a discrete analog by considering a 
quadrature formula in place of the integral. 

The discretization process requires to change the curve y(t) in the phase space R 2m 
to a simpler curve a(t) (generally but not necessarily a polynomial), which is meant 
to yield the approximation at time t\ = to + h, that is y(to + h) = a(to + h) + 0(h p+1 ), 
where p is the order of the resulting numerical method. In a certain sense, the problem 
of numerically solving (1.1) while preserving the Hamiltonian function is translated 
into a quadrature problem. 

For example, consider the segment a(to + ch) = (1 — c)yo + cyi, with c <G [0, 1], 
joining yo to an unknown point y\ of the phase space. The line integral of V-ff(y) 
evaluated along a becomes 

H( yi ) - H(yo) = %j - y f f VfT((l - c)y + c Vl )dc. (1.3) 

Jo 

Now assume that H(y) = H(q,p) is a polynomial of degree v in the generalized 
coordinates q and in the momenta p. The integrand in (1.3) is a polynomial of 
degree v — 1 in c and can be exactly solved by any quadrature formula with abscissae 
ci < ca < • • • < Ck in [0, 1] and weights b\, . . . , bk, having degree of precision d > v — \. 
We thus obtain 

k 

H{ Vl ) - H(y ) = h{yx - yofY. b ^ H (( l - c i)Vo + cm)- 

i=l 

To get the energy conservation property we impose that y\ — yo be orthogonal to the 
above sum, and in particular we choose (for the sake of generality we use f(y) in place 
of J"VH(y) to mean that the resulting method also makes sense when applied to a 
general ordinary differential equation y' = f(y)) 

k 

yi = yo + h'y y j bif(Y i ), Yl = (1 - cijyo + Ciy-L, i=l,...,k. (1.4) 



Formula (1.4) defines a Runge-Kutta method with Butcher tableau — 
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-pi — , where 

c and b are the vectors of the abscissae and weights, respectively. The stages Yi are 
called silent stages since their presence does not affect the degree of nonlinearity of 
the system to be solved at each step of the integration procedure: the only unknown 
is yi and consequently (1.4) defines a mono-implicit method. Mono-implicit methods 
of Runge-Kutta type have been researched in the past by several authors (see, for 
example, [9, 1, 10, 8] for their use in the solution of initial value problems). 

Methods such as (1.4) date back to 2007 [18, 19] and are called fc-stage trapezoidal 
methods since on the one hand the choice k = 2, c\ ~ 0, C2 = 1 leads to the trapezoidal 
method and on the other hand all other methods evidently become the trapezoidal 
method when applied to linear problems. 



Generalizations of (1.4) to higher orders require the use of a polynomial a of 
higher degree and are based upon the same reasoning as the one discussed above. 
Up to now, such extensions have taken the form of Rungc-Kutta methods [4, 5, 
6]. It has been shown that choosing a proper polynomial a of degree s yields a 
Runge-Kutta method of order 2s with k > s stages. The peculiarity of such energy- 
preserving formulae, called Hamiltonian Boundary Value Methods (HBVMs), is that 
the associated Butcher matrix has rank s rather than k, since k — s stages may be 
cast as linear combinations of the remaining ones, similarly to the stages Yi in (1.4). 1 
As a consequence, the nonlinear system to be solved at each step has dimension 2ms 
instead of 2mk, which is better visualized by recasting the method in block-BVM 
form [4]. 

In the case where H(y) is not a polynomial, one can still get a practical energy 
conservation by choosing k large enough so that the quadrature formula approximates 
the corresponding integral to within machine precision. Strictly speaking, taking the 
limit as k — > oo leads to limit formulae where the integrals come back into play in 
place of the sums. For example, letting k — > oo in (1.4) just means that the integral 
in (1.3) must not be discretized at all, which would yield the Averaged Vector Field 
method y\ = yo + h f Q /((l — c)yo + cy\) dc, (see [11, 24] for details). 

In this paper we start an investigation that follows a different route. Unlike 
the case with HBVMs, we want now to take advantage of the previously computed 
approximations to extend the class (1.4) in such a way to increase the order of the 
resulting methods, much as the class of linear multistep method may be viewed as a 
generalization of (linear) one step methods. The general question we want to address 
is whether there exist fc-step mono- implicit energy-preserving methods of order greater 
than two. Clearly, the main motivation is to reduce the computational cost associated 
with the implementation of HBVMs. 

The purpose of the present paper is to give an affermative answer to this issue in 
the case k = 2. More specifically, the method resulting from our analysis, summarized 
by formula (4.1), may be thought of as a nearly linear two-step method in that it is 
the sum of a fourth order linear two-step method, formula (4.3), plus a nonlinear 
correction of higher order. 

The paper is organized as follows. In Section 2 we introduce the general formula- 
tion of the method, by which we mean that the integrals are initially not discretized 
to maintain the theory at a general level. In this section we also report a brief de- 
scription of the HBVM of order four, since its properties will be later exploited to 
deduce the order of the new method: this will be the subject of Section 3. Section 4 is 
devoted to the discretization of the integrals, which will produce the final form of the 
methods making them ready for implementation. A few test problems are presented 
in Section 5 to confirm the theoretical results. 

2. Definition of the method. Suppose that y\ is an approximation to the true 
solution y(t) at time tx = to + h, where h > is the stepsize of integration. More 
precisely, we assume that 

(At) y(ti) = J/i + 0{hP +1 ) with p > 4; 

(A2) H(yi) = H(yo), which means that y\ lies on the very same manifold H(y) = 
H(yo) as the continuous solution y(t). 



A documentation about HBVMs, Matlab codes, and a complete set of references is available at 
the url [3]. 
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The two above assumptions arc fulfilled if, for example, we compute y\ by means of 
a HBVM (or an oo-HBVM [5]) of order p > 4. The new approximation y 2 — y{t%) = 
y(to + 2h) is constructed as follows. 

Consider the quadratic polynomial a (to + 2rh) that interpolates the set of data 
{(to +jh, yj)}j=o,i.2- Expanded along the Newton basis {Pj(r)} defined on the nodes 
To = 0, T\ = h, T2 = 1, the polynomial a takes the form (for convenience we order the 
nodes as To , t 2 , ri ) 

a(to + 2rh) = y + (y 2 - Vo)r + 2{y 2 - 2 Vl + y )r(r - 1). (2.1) 

As r ranges in the interval [0, 1], the 2m-length vector 7(7-) = a(to + 2rh) describes a 
curve in the phase space R 2m . The line integral of the conservative vector field VH(y) 
along the curve 7 will match the variation of the energy function H(y), that is 



H(y 2 )-H(y ) = f WH(y)dy= f [i(r)] T Vif( 7 (r)) dr 



{V2 - yo)' 



' [ WH( 1 (T))dT + 2(y 2 -2y 1 +y Q ) T f (2r - 1)VH( T (t)) dr. 
Jo Jq 



The energy conservation condition H(y 2 ) = H(yo) yields the following equation in 
the unknown z = y 2 



(z - yo) : 



/ VH( 7 (r))dr = -2(z - 2 Vl + y ) T [ (2r - 1)V# ( 7 (r)) dr. (2.2) 
Jo ^0 



The method we are interested in has the form y 2 — ^hijJOi j/i ) , where is implicitly 
defined by the following nonlinear equation in the unknown z: 



z = y + 2hJa(z) + u JJ u2 a(z), with a(z) = / Vff( 7 (r)) dr, (2.3) 



\\°>(z)\\2 

where the residual r(z) is defined as 



f Vff( 7 (r))dr, 
Jo 



r{z) = -2(z - 2 Vl + yo) 1 ' [ (2r - l)Vi7( 7 (r)) dr. (2.4) 

Jo 

A direct computation shows that any solution z* of (2.3) also satisfies (2.2). In the 
next section we will show that (2.3) admits a unique solution y 2 = z* satisfying the 
order condition y 2 = y(to + 2h) + 0(h 5 ). Such a result will be derived by regarding 
(2.3) as a perturbation of the HBVM of order 4 and, in turn, by comparing the 
two associated numerical solutions. To this end and to better explain the genesis of 
formula (2.3) and the role of the integrals therein, a brief introduction of the HBVM 
formula of order four is in order. 

2.1. HBVM of order four. Suppose that both j/i and y 2 are unknown (so now 
yi is no longer given a priori as indicated by assumption (A\)): let us call them ui 
and u 2 respectively. For (2.2) to be satisfied, we can impose the two orthogonality 
conditions 



U2-yo = VihJ / V-ff(7(r)) dr, 

u 2 - 2u x +yo= mhJ [ (2r - l)Vi?( 7 (r)) dr, 
Jo 



(2.5) 



giving rise to a system of two block-equations (the curve 7(7") = o~(to + 2rh) is as in 
(2.1) with U\ and u 2 in place of y\ and y 2 )- Setting the free constants r\\ and 772 equal 
to 2 and 3, respectively, confers the highest possible order, namely 4, on the resulting 
method: u 2 = y(to + 2h) + 0(h 5 ) (see [19] for details). 2 Furthermore, it may be shown 
that the internal stage u\ satisfies the order condition u\ = y(to + h) + 0(h A ). 

Evidently, the implementation of (2.5) on a computer cannot leave out of consid- 
eration the issue of solving the integrals appearing in both equations. Two different 
situations may emerge: 

(a) the Hamiltonian function H (y) is a polynomial of degree v. In such a case, the 
two integrals in (2.5) are exactly computed by a quadrature formula having 
degree of precision d > 2v — 1. 

(b) H(y) is not a polynomial, nor do the two integrands admit a primitive func- 
tion in closed form. Again, an appropriate quadrature formula can be used 
to approximate the two integrals to within machine precision, so that no sub- 
stantial difference is expected during the implementation process by replacing 
the integrals by their discrete counterparts. 

Case (a) gives rise to an infinite family of Runge-Kutta methods, each depending 
on the specific choice (number and distribution) of nodes the quadrature formula is 
based upon (see [5] for a general introduction on HBVMs and [6] for their relation with 
standard collocation methods). For example, choosing k nodes according to a Gauss 
distribution over the interval [0, 1] results in a method that precisely conserves the 
energy if applied to polynomial canonical Hamiltonian systems with v < k and that 
becomes the classical 2-stage Gauss collocation method when k = 2. On the other 
hand, choosing a Lobatto distribution yields a Runge-Kutta method that preserves 
polynomial Hamiltonian functions of degree v < k — 1 and that becomes the Lobatto 
III A method of order four when k = 2. 

The method resulting from case (b) are undistinguishable from the original for- 
mulae (2.5) in that they are energy-preserving up to machine precision when applied 
to any regular canonical Hamiltonian system. Stated differently, (2.5) may be viewed 
as the limit of the family of HBVMs of order four, as the number of nodes tends to 
infinity. For this reason the limit formulae (2.5) have been called oo-HBVM of order 
4 (see [5]). 

Remark 1. In the present context, y\ being a known quantity, the unknown z 
in (2.2) cannot in general satisfy, at the same time, both orthogonality conditions in 
(2.5). However, since y\ may be thought of as an approximation of order four to the 
quantity u\ in (2.5), should we only impose the first orthogonality condition, namely 

z -y =2hJa(z), (2.6) 

we would expect the residual r(z) (the right hand side of (2.2) ) to be very small. 3 This 
suggests that a solution to (2.2) that yields an approximation of high order to y(to+2h) 
may be obtained by allowing a small deviation from orthogonality in (2.6). This is 
accomplished by setting z — j/o = 2hJa{z) + Sa(z), and by tuning the perturbation 
parameter 5 in such a way that (2.2) be satisfied: this evidently gives 5 = | | | ^ and 
we arrive at (2.3). 



2 Since we are integrating the problem on an interval [to,t2] of length 2h, we have scaled the 
constants vi and r]2 by a factor two with respect to the values reported in [19], 

3 By exploiting the result in Lemma 3.2 below, it is not difficult to show that actually (2.6) implies 
r(z) = 0(h 5 ). This aspect is further emphasized in the numerical test section (sec Table 5.2). 
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3. Analysis of the method. Results on the existence and uniqueness of a 
solution of (2.3) as well as on its order of accuracy will be derived by first analyzing 
the simpler nonlinear system 



yo + 2hJa(z), with a(z) 



f VfT(7(T))dT, (3.1) 
Jo 



obtained by neglecting the correction term -^^ I a{z). For z 6 K 2m we set (see (2.1)) 
1z{t) =y + {z- y )T + 2(z - 2yi + y )r(r - 1), (3.2) 

and (see (3.1)) 

$(z) = y + 2hJa(z). (3.3) 

In the following 1 1 • | will denote the 2-norm. 

Lemma 3.1. There exist positive constants p and ho such that, for h < ho, system 
(3.1) admits a unique solution z in the ball B(yo,p) of center yo and radius p. 

Proof. We show that constants h ,p > exist such that the function defined in 
(3.3) satisfies the following two conditions for h < ho: 

(a) $(z) is a contraction on B(yo,p), namely 

Vz,w 6 B(y ,p), \\$(z) — ®(w)\\ < L\\z — w\\, with L < 1; 

(b) ||$(y )-yo|| < (l-L)p. 

The contraction mapping theorem can then be applied to obtain the assertion. 

Let B(yo,p) a ball centered at yo with radius p. We can choose h' and p small 
enough that the image set il = {7z(t) : r <E [0, 1], z £ B(yo,p), h < h' } is entirely 
contained in a ball B(yo,p') which, in turn, is contained in the domain of V 2 iJ(y). 4 
We set 

M p = max ||V 2 ffH||. 
From (3.1) and (3.2) we have 

^ = J V 2 H{ 1z {t))^At = J V 2 if( 7 ,(r)) r(2r - 1) dr 

and hence 

da(z) 



dz 

Consequently (a) is satisfied by choosing 



< M p I t\2t- l|dr = -M n . 



L = \M P (3.4) 



and ho < minj-p-, h' }- Concerning (b), we observe that 

*(l/o) - Vo = 2hJa(y ) = 2hJ f WH(y + 4(y - Vi )t{t - 1)) dr, 



Notice that, by definition, the set f2 is an open simply connected subset of R 2m containing 
B(yo,p) while, from the assumption (Ai), decreasing h causes the point yi to approach yo. 

6 



hence ||$(yo 



Vol 



2h\\a(y )\\ with ||a(yo)|| bounded with respect to h. Since L 



vanishes with h (see (3.4)), we can always tune ho in such a way that 2/i||a(y )| < 
(1 - L)p. □ 

Lemma 3.2. The solution z of (3.1) satisfies y(t + 2h) — z = 0(h 5 ). 

Proof. Under the assumption (Ai), (3.1) may be regarded as a perturbation of 
system (2.5), since y\ and u\ are 0(h 5 ) and 0(h A ) close to y(t + h) respectively. 5 
Since u 2 = y(t + 2h) + 0(h 5 ), we can estimate the accuracy of z as an approximation 
of y(t + 2h) by evaluating its distance from u-x- 

Let j(t) be the underlying quadratic curve associated with the HBVM defined 
by (2.5), namely 

7(t) = yo + (u 2 - yo)r + 2(u 2 - 2u x + Uq)t(t - 1). (3.5) 
Considering that (see (3.2)) 

lu 2 (r) = y + (u 2 - y )r + 2(u 2 - 2y t + yo)r(T - 1) = j(t) + 4(u x - yi)r(r - 1), 
from the first equation in (2.5) and (3.3) we get 



$(u 2 ) 



yo + 2hJ f VH{ lu2 (r)) dr = y + 2hJ f VH(j(t)) dr 
Jo Jo 

+ShJ f V 2 H(j(t))t(t -l)dT-(«i-yi) + 0(||ui-|/i|| a ) 
Jo 

= u 2 + 0(h 5 ). 

If h is small enough, u 2 will be inside the ball B(yo,p) defined in Lemma 3.1. The 
Lipschitz condition yields (see (3.4)) 

\\z-u 2 \\ = ||$(z) - $(u 2 ) + 0{h 5 )\\ < -M p \\z ~ u 2 \\ + 0(/i 5 ), 

and hence \\z — u 2 \ \ = 0(h 5 )\\. □ 

The above result states that (3.1) defines a method of order 4 which is a simplified 
(non corrected) version of our conservative method defined at (2.3). In Section 5 the 
behavior of these two methods will be compared on a set of test problems. We now 
state the analogous results for system (2.3). 

Theorem 3.3. Under the assumption (A\), for h small enough, equation (2.3) 
admits a unique solution z* satisfying y(t + 2h) — z* = Olh 5 ). 

Proof. Consider the solution z of system (3.1). We have (see (3.5)) 



7«W - = {z- u 2 )t{2t - 1) + 4(«i - yi)r(r - 1) = 0{h 5 ), 



and 



z - 2yi + y =u 2 - 2ui + yo + 0(h 5 ). 
Hence, by virtue of (2.5), 

i T 



r(z) = -2[(u 2 -2 Ul +y a ) + 0(h 5 )]- 



(2t - 1)V#(7(t)) dr + 0(h 5 ) 



0(h 5 ). 



5 This also implies that u\ — y\ = 0(h 4 ). 



Since a{z) is bounded with respect to h, it follows that, in a neighborhood of z, system 
(2.3) may be regarded as a perturbation of system (3.1), the perturbation term being 

Consider the ball B(z, R(z, h)): since z = y + 0(h), and R(z,h) = 0(h 5 ), 
this ball is contained in B(yo,p) defined in Lemma 3.1 and the perturbed function 
$(z) + R(z,h) is a contraction therein, provided h is small enough. Evaluating the 
right-hand side of (2.3) at z = z we get 

y + 2hJa(z) + R(z, h) = z + R(z, h), 

which means that property (b) listed in the proof of Lemma 3.1, with z in place of 
yo, holds true for the perturbed function yo + 2hJa(z) + R(z, h), and the contraction 
mapping theorem may be again exploited to deduce the assertion. □ 

4. Discretization. As was stressed in Section 2, formula (2.3) is not operative 
unless a technique to solve the two integrals is taken into account. The most obvious 
choice is to compute the integrals by means of a suitable quadrature formula which 
may be assumed exact in the case where the Hamiltonian function is a polynomial, 
and to provide an approximation to within machine precision in all other cases. 

Hereafter we assume that H(q,p) is a polynomial in q and p of degree v. Since 
j(t) has degree two, it follows that the integrand functions appearing in the definitions 
of a(z) and r(z) at (2.3) and (2.4) have degree 2v — 2 and 2v — \ respectively and can 
be solved by any quadrature formula with abscissae c\ < c% < • • • < c& in [0, 1] and 
weights b\, . . . , bk, having degree of precision d > 2u — 1. In place of (2.3) we now 
consider the equivalent form suitable for implementation 

k 

y 2 = yo + 2hjJ2b l VH(j{c l )) + G(y , yi ,y 2 ), (4.1) 

i=l 

where 

rl , -2(y 2 -2y 1 +y ) T j:t 1 b i (2c l -l)VH( 1 (c l )) ^ UV7ZJ , , 
G(y ,yi,y2) = k ^b^H^)). 

IIEi=i&iVfr(7(ci))|^ i=1 

Notice that from (2.1) we get 

7(c<) = (1 - 3c H + 2c})y Q + 4c,(l - a) yi + Cl (2c, - l)y 2 , (4.2) 

that is, "f(ci) is a linear combination, actually a weighted average, of the approxi- 
mations yo, yi and y 2 . Therefore, since G(yo, j/i, y 2 ) = 0(h b ) (see Lemma 3.2 and 
Theorem 3.3), we may look at this term as a nonlinear correction of the generalized 
linear multistep method 

fc 

y 2 = y a + 2hjJ2kVH( 7 (ci)). (4.3) 

i=i 

Example 1. If H(q,p) is quadratic, we can choose k = 3, c\ = 0, c% = ^> C3 = 1, 
bi = &3 = i and b 2 = |, that is we can use Simpson's quadrature formula to compute 
the integrals in (2.3) and (2.4). Since, in such a case, 7(0,) = yi—i, method (4.3) 
becomes 



V2 = yo + 3 J (Vi?(2/ ) + 4V£% 1 ) + VH{y 2 )) , 
8 



that is, the standard Milne-Simpson's method. 

In all other cases 7(0$) will differ in general from yj, j = 1,2,3 and may be 
regarded as an off-point entry in formula (4.3). In the sequel we will denote the 
method defined at (4.1) by Mj. and its linear part, defined at (4.3), by M' k . Of course, 
the choice of the abscissae distribution influences the energy preserving properties of 
the method Mk, as is indicated in Table 4.1. 



Abscissae distribution: uniform 


Lobatto 


Gauss 


Energy preserving when: degH < |~|] 


degH <k—l 


deg H < k 



Table 4.1 



Energy preserving properties of method Mfc for some well-known distributions of the nodes {ci}- 



5. Numerical tests. Hereafter we implement the order four method Mk on a 
few Hamiltonian problems to show that the numerical results are consistent with the 
theory presented in Section 3. In particular, in the first two problems the Hamiltonian 
function is a polynomial of degree three and six respectively, while the last numerical 
test reports the behavior of the method on a non-polynomial problem. 

Each step of the integration procedure requires the solution of a nonlinear system, 
in the unknown j/2, represented by (4.1) for the method Mk and (4.3) for the method 
M' k . The easiest way (although not the most efficient one) to find out a solution is by 
means of fixed point iteration that, in the case of the method Mk, reads 

k 

z s+1 = y + 2hJj2 b iVH('y Zs (c i )) + G(yo,yuz s ), 8 = 1,2,..., (5.1) 

i=l 

where 7 2 is defined at (3.2) and Zq is an initial approximation of yi which is then 
refined by setting ?/2 = z s with z s ~ lim s _ ! . 00 z s . From Theorem 3.3 and the preceding 
lemmas we deduce that such a limit always exists provided that h is small enough. The 
value of zq could be retrieved via an extrapolation based on the previous computed 
points or by considering the method M' k as a predictor for Mk- 

We will consider a Lobatto distribution with an odd number k of abscissae {ct}. 
In fact, if k is odd, since yo = 7(0) = j(c\) and y\ = 7(5) = 7(C|-*-|), we save two 
function evaluations during the iteration (5.1). 

5.1. Test problem 1. The Hamiltonian function 

H(q,p) = \p 2 + \q 2 - jU 3 (5.2) 

defines the cubic pendulum equation. We can solve it by using five Lobatto nodes 
to discretize the integrals in (2.3), thus getting the method M5. The corresponding 
numerical solution, denoted by (q n ,Pn), is plotted in Figure 5.1. For comparison 
purposes we also compute the numerical solution (q' n ,p' n ) provided by the fourth 
order method, say M5, obtained by neglecting in (2.3) the correction term, that is 
by posing r(z) = 0. Figure 5.2 clearly shows the energy conservation property, while 
Table 5.1 summarizes the convergence properties of the two methods. 

5.2. Test problem 2. The Hamiltonian function 

H(p, q) = -p 3 --p+ —q 6 + -q 4 - -q 3 + - (5.3) 
9 



100 200 300 400 500 600 -1 -05 0.5 1 1.5 

Fig. 5.1. Numerical solution (q n ,Pn) versus time t n (left picture) and on the phase plane (right 
picture). Parameters: initial condition yo = [0, 1]; stepsize h = 0.5; integration interval [0,2007r], 
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Fig. 5.2. Hamiltonian function evaluated along the numerical solution (pn,9n) (horizontal line) 
and along the numerical solution (p' n ,q' n ) (irregularly oscillating line). 

has been proposed in [14] to show that symmetric methods may suffer from the energy 
drift phenomenon even when applied to reversible systems, that is when H(—p, q) = 
H(p,q). e For our experiment, we will use i/q = [0.2, 0.5] as initial condition. 

Since dcg(H(q,p)) = 6, we need a Lobatto quadrature based on at least seven 
nodes to assure that the integrals in (2.3) are computed exactly. Therefore we solve 
(5.3) by method My. For comparison purposes, it is also interesting to show the 
dynamics of the symmetric non-conservative method M 7 . Figure 5.3 displays the 
results obtained by the two methods implemented with stepsize h = over the 
interval [0, 10 3 ]. In particular, the numerical trajectories generated by method Af 7 
and M7, are reported in the left-top and left-bottom pictures respectively, while the 
right picture reports the corresponding error in the Hamiltonian function evaluated 
along the two numerical solutions, namely \H{y n ) — H(yo)\. 

Evidently, the numerical solution produced by M 7 rapidly departs from the level 
curve H(q,p) — H(qo,Po) but it remain eventually bounded and the points (q n ,Pn) 
seem to densely fill a bounded region of the phase plane. 

On the contrary, since the degree of freedom of the present problem is one, the 
points (q n ,p n ) produced by M7 lie on the very same continuous trajectory covered by 



6 In fact, the authors show that the system deriving from (5.3) is equivalent to a reversible system 
(see also [7, 22] for a discussion on the integration of reversible Hamiltonian systems by symmetric 
methods). 
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Table 5.1 

Methods M5 (with correction term) and Mg (without correction term) are implemented on 
the cubic pendulum equation (5.2) on the time interval [0, 10] {or several values of the stepsize h. 

error ( — ) 

The order of convergence is numerically evaluated by means of the formula log 2 crror ^ ■ As was 
expected, the maximum displacement of the numerical Hamiltonian H(y n ) from the theoretical value 
H{yo) is close to the machine precision for the method M5, independently of the stepsize h used. 



y(t): this is also confirmed by looking at the bottom graph in the right picture. 

Table 5.2 shows the behavior of method M7 applied to problem (5.3) as the 
stepsizes h goes to zero. Notice the 0(h 5 ) rate of convergence to zero for the residual 
function r(z) in (2.4). 
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Fig. 5.3. Left pictures: numerical solutions in the phase plane computed by method Afy (top 
picture) and M7 (bottom picture). Right picture: error in the numerical Hamiltonian function 
\H{Vn) — H(yo)\ produced by the two methods. Parameters: initial condition yo = [0.2, 0.5]; stepsize 
h = 0.1; integration interval [0, 1000]. 



5.3. Test problem 3. We finally consider the non-polynomial Hamiltonian 
function 

H(q 1 ,q 2 ,Pi,p 2 ) = +Pl) ~ , } a (5.4) 
L V <7i + <7 2 

that defines the well known Kepler problem, namely the motion of two masses under 
the action of their mutual gravitational attraction. Taking as initial condition 

(gi(0), 92 (0),pi(0),p 2 (0))= ( l-e, 0, 0, (5.5) 
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Table 5.2 

Performance of method M-j applied to problem (5.3), with initial condition yo = [0.2, 0.5], on 
the time interval [0, 250] for several values of the stepsize h, as specified in the first column. The 
second and third columns report the relative error in the last computed point j/jv, N = T/h and the 
corresponding order of convergence. Since the integrals appearing in (2.3) are precisely computed by 
the Lobatto quadrature formula with seven nodes, the error in the numerical Hamiltonian H{yj^) is 
zero up to machine precision. The last two columns list the residual r(yjv) defined in (2.4) and its 
order of convergence to zero. 



yields an elliptic periodic orbit of period 2tt and eccentricity e G [0,1). We have chosen 
e = 0.6. Though the vector field fails to be a polynomial in q\ and q 2 , we can plan 
to use a sufficiently large number of quadrature nodes to discretize the integrals in 
(2.3) so that the corresponding accuracy is within the machine precision. Under this 
assumption, and taking aside the effect of the floating point arithmetic, the computer 
will make no difference between the conservative formulae (2.3) and their discrete 
counterparts. 

The left picture in Figure 5.4 explains the above argument. It reports the error 
\H(y n ) — H(yo)\ in the Hamiltonian function for various choices of the number of 
Lobatto nodes, and precisely k = 3, 5, 7, 9. We see that the error decreases quickly 
as the number of nodes is incremented and for k = 9 it is within the epsilon machine. 7 

The use of finite arithmetic may sometimes cause a mild numerical drift of the 
energy over long times, like the one shown in the upper line in the right picture 
of Figure 5.4. This is due to the fact that on a computer the numerical solution 
satisfy the conservation relation H(y n ) = H(yo) up to machine precision times the 
conditioning number of the nonlinear system that is to be solved at each step. 

To prevent the accumulation of roundoff errors we may apply a simple and costless 
correction technique on the approximation y n which consists in a single step of a 
gradient descent method (see also [2]). More precisely, the corrected solution y* n is 
defined by 

Vg(yn) ... H{y n ) - H(y ) 

y n = yn-a—— — -—, with a = — - — , 5.6 

\\VH(y n )\\ 2 \\VH(y n )\\ 2 

which stems from choosing as a the value that minimizes the linear part of the function 
F(a) — H(y n — a ^^^"^ ) — H(yo). The bottom line in the right picture of Figure 
5.4 shows the energy conservation property of the corrected solution. 

6. Conclusions. We have derived a family of mono-implicit methods of order 
four with energy-preserving properties. Each element in the family originates from a 



7 A11 tests were performed in Matlab using double precision arithmetic. 
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Fig. 5.4. Left picture. Error in the numerical Hamiltonian function \H(y„) — H(yo)\ produced 
by methods M^, with k = 3, 5, 7, 9. Parameters: stepsize h = 0.05, integration interval [0, 50] . Right 
picture. Roundoff errors may cause a drift of the numerical Hamiltonian function (upper line) which 
can be easily taken under control by coupling the method with a costless correction procedure like the 
one described at (5.6). 



limit formula and is denned by discretizing the integral therein by means of a suitable 
quadrature scheme. This process assures an exact energy conservation in the case 
where the Hamiltonian function is a polynomial, or a conservation to within machine 
precision in all other cases, as is also illustrated in the numerical tests. Interestingly, 
each method may be conceived as a 0(h 5 ) perturbation of a two-step linear method. 
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